Differential Equations

Adapted from Woflram.com documentation

Automatically selecting between hundreds of powerful and in many cases original algorithms, the Wolfram Language provides both numerical and symbolic solving of differential equations (ODEs, PDEs, DAEs, DDEs, ...). With equations conveniently specified symbolically, the Wolfram Language uses both its rich set of special functions and its unique symbolic interpolating functions to represent solutions in forms that can immediately be manipulated or visualized

Download original notebook

Take derivatives

Take a derivative symbolically of any function

Sin'[x]

Cos[x]

This is an equvivalent to

D[Sin[x], x]

Cos[x]

Define your own function in a separate symbol

f[x_] := Sin[x] + (*SpB[*)Power[x(*|*),(*|*)2](*]SpB*)
f'[x]
f'[0.5]

2 x+Cos[x]

1.8775825618903728`

DSolve

Solve a differential equation with dependent variable y[x]:

DSolve[y'[x] + y[x] == a Sin[x], y[x], x]

{{y[x]->((*SpB[*)Power[E(*|*),(*|*)-x](*]SpB*)) (*SbB[*)Subscript[C(*|*),(*|*)1](*]SbB*)+((*FB[*)((1)(*,*)/(*,*)(2))(*]FB*)) a (-Cos[x]+Sin[x])}}

Include a boundary condition:

DSolve[{y'[x] + y[x] == a Sin[x], y[0] == 0}, y[x], x]

{{y[x]->-(*FB[*)((1)(*,*)/(*,*)(2))(*]FB*) a ((*SpB[*)Power[E(*|*),(*|*)-x](*]SpB*)) (-1+((*SpB[*)Power[E(*|*),(*|*)x](*]SpB*)) Cos[x]-((*SpB[*)Power[E(*|*),(*|*)x](*]SpB*)) Sin[x])}}

Specify only the first argument of DSolve:

{{y[x]->((*SpB[*)Power[E(*|*),(*|*)x](*]SpB*)) (*SbB[*)Subscript[C(*|*),(*|*)1](*]SbB*)}}

ODE

Solve Riccati equaiton

DSolve[{y'[x] == y[x] (1 - y[x]/27), y[0] == a}, y, x] // Quiet

{{y->Function[{x},(*FB[*)((27 a ((*SpB[*)Power[E(*|*),(*|*)x](*]SpB*)))(*,*)/(*,*)(27-a+a ((*SpB[*)Power[E(*|*),(*|*)x](*]SpB*))))(*]FB*)]}}

Plot the solution for different initial values:

Plot[Evaluate[y[x] /. % /. {{a->(*FB[*)((1)(*,*)/(*,*)(13))(*]FB*)},{a->(*FB[*)((1)(*,*)/(*,*)(2))(*]FB*)},{a->4}}], {x, 
  0, 18}]

(*VB[*)(FrontEndRef["4cefbc52-2c0e-4587-8f4d-245e34232e5c"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKmySnpiUlmxrpGiUbpOqamFqY61qkmaToGpmYphqbGBkbpZomAwCLthWk"*)(*]VB*)

Solve a linear pendulum equation:

sol = DSolve[{y''[x] + y[x] + (*FB[*)((1)(*,*)/(*,*)(5))(*]FB*) y'[x] == 0, y[0] == 1, y'[0] == 1/3}, y , x]

{{y->Function[{x},((*FB[*)((1)(*,*)/(*,*)(99))(*]FB*)) ((*SpB[*)Power[E(*|*),(*|*)-x/10](*]SpB*)) (99 Cos[((*FB[*)((3)(*,*)/(*,*)(10))(*]FB*)) ((*SqB[*)Sqrt[11](*]SqB*)) x]+13 ((*SqB[*)Sqrt[11](*]SqB*)) Sin[((*FB[*)((3)(*,*)/(*,*)(10))(*]FB*)) ((*SqB[*)Sqrt[11](*]SqB*)) x])]}}

Plot[y[x] /. sol, {x, 0, 27}]

(*VB[*)(FrontEndRef["81c57fed-b6c9-4349-9c76-830e52d5e92d"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKWxgmm5qnpaboJpklW+qaGJtY6lomm5vpWhgbpJoapZimWhqlAACFNxWP"*)(*]VB*)

Animate it

{x, y[If[x < t, x, t]]} /. First[sol];
AnimateParametricPlot[%, {x, 0, 40}, {t,0, 40}, PlotRange->{{0,40}, {-1,1}}]

(*VB[*)(FrontEndRef["76dcec58-904e-4952-9d94-7bb0b20a968b"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKm5ulJKcmm1roWhqYpOqaWJoa6VqmWJromiclGSQZGSRamlkkAQCDvxWF"*)(*]VB*)

Model a block on a moving conveyor belt anchored to a wall by a spring. Compare positions and velocities for the different values of the parameters of the system (mass of the block, belt speed, friction coefficient, spring constant):

ClearAll[m];

beltv[t_] = vb;
spring[x_] = k (l - x);
friction[v_] := -a (v - beltv[t]);

Newton's equation for the block:

sys := {
  m x''[t] == spring[x[t]] + friction[x'[t]], 
  x[0] == l, 
  x'[0] == 0
};

Solve for position and velocity:

sol = DSolve[sys, x, t];

{pos[m_, k_, a_, vb_, t_], vel[m_, k_, a_, vb_, t_]} = {x[t], x'[t]} /. sol[[1]] // FullSimplify;

The block stabilizes just above the spring's natural length of l:

l = 1;
Manipulate[{Plot[{pos[m, k, a, vb, t]}, {t, 0, 2}, PlotRange -> All, 
   PlotLabel -> "Position"], 
  Plot[{vel[m, k, a, vb, t], vb}, {t, 0, 2}, PlotRange -> All, 
   PlotLabel -> "Velocity"]} // Row, {{m, 4, "mass of the block"}, 4, 10, 
  4}, {{k, 250, "spring constant"}, 250, 1000, 
  250}, {{vb, 5, "belt velocity"}, 5, 10, 
  5}, {{a, 15, "friction coefficient"}, 5, 35, 15}]

(*GB[*){{(*VB[*)(EventObject[<|"Id" -> "7172d3a7-25ab-4438-ac8f-1cdbea0c19c1", "Initial" -> {4, 250, 5, 15}, "View" -> "abe8f632-7809-45ec-bbe3-b1d047b75f45"|>])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKJyalWqSZGRvpmlsYWOqamKYm6yYlpRrrJhmmGJiYJ5mbppmYAgCLnhXY"*)(*]VB*)}(*||*),(*||*){(*VB[*)(FrontEndRef["61febc24-f2fb-4ebc-94e2-7421f388c458"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKmxmmpSYlG5nophmlJemaANm6liapRrrmJkaGacYWFskmphYAkjcV3Q=="*)(*]VB*)}}(*]GB*)

Hybrid Differential Equations

Model a damped oscillator that gets a kick at regular time intervals:

system = {x''[t] + .1 x'[t] + x[t] == 0, x[0] == 1, x'[0] == 0};
control = WhenEvent[Mod[t, 1] == 0, x'[t] -> x'[t] + 1];

sol = DSolve[{system, control}, x, {t, 0, 50}];
Plot[Evaluate[{x[t], x'[t]} /. sol], {t, 0, 50}, PlotLegends->LineLegend[Automatic, {"Posiiton", "Velocity"}]]

(*VB[*)(Legended[ToExpression[FrontEndRef["545f6ca9-eab4-45e6-8832-cedc643764b2"], InputForm], Placed[LineLegend[{Directive[Opacity[1.], RGBColor[0.24, 0.6, 0.8], AbsoluteThickness[2]], Directive[Opacity[1.], RGBColor[0.95, 0.627, 0.1425], AbsoluteThickness[2]]}, {"Posiiton", "Velocity"}, LegendMarkers -> None, LabelStyle -> {}, LegendLayout -> "Column"], After, Identity]])(*,*)(*"1:eJylUUtOwzAQDX+oQAIOgITENps2DWVRRYUWhJTS0lTdO8kYrAa7sh0gB+AQsGXHCbrmACzYcADYICQEN8B2+ChsWDCLJ3v8Zub5zXrIenjKsiyxpGBA4KwJEeNIMh7MqYwPR0DjMp7UlJKCVkzUmybiCZ1bVbDLGZUtGrfOIUolChMINlS66lSxG6EtG1Do2E4VXLtWq5TtCOLIdSqbrhN+Np5W0EtV2bw+AIo7NMlMts9TyPXNKugmSBXjmS8xPqGQK/zp4xMh84oFBU3CIZLkFHK1+kudEYqIzLhl4s3LyWby3vYOSxjn47WLl8PxnccrJh49fnWp49nL26woaISCJamE/jGJhhSEIFrCfydjE68ev3m/b4fLTx6vlx6uR/XbvycXHDBOdpkgalvUXAaQMD296LhZe25hG/EhcGGeDhiFX0TjNwohCWSWALYKfhepi989fZSxVAZ6d+p/6QkVencNLNUgLWo/BiqVqA+5lqbU"*)(*]VB*)

In a square box, model a ball that changes direction upon impact with the side walls:

ClearAll[a,b,x]

sol = DSolve[{x'[t] == a[t], y'[t] == b[t], x[0] == 0, y[0] == 0, 
    a[0] == 1, b[0] == 17/12, WhenEvent[x[t]^2 == 1, a[t] -> -a[t]], 
    WhenEvent[y[t]^2 == 1, b[t] -> -b[t]]}, {x, y}, {t, 0, 50}, 
   DiscreteVariables -> {a, b}];

ParametricPlot[{x[t], y[t]} /. sol, {t, 0, 50}, Frame -> True, 
 FrameTicks -> None, Axes -> False]

(*VB[*)(FrontEndRef["c14b8965-b616-4a8d-8d2e-49751f590327"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKJxuaJFlYmpnqJpkZmumaJFqk6FqkGKXqmliamxqmmVoaGBuZAwB86BT2"*)(*]VB*)

NDSolve

Solve a first-order ordinary differential equation:

s = NDSolve[{y'[x] == y[x] Cos[x + y[x]], y[0] == 1}, y, {x, 0, 30}]

{{y->(*VB[*)BoxForm`temporalStorage$506528(*,*)(*"1:eJxdkUeP2lAAhEmTovyKZE9gIrnwMBBpD/biXhbjgp3VHmz8MO7tGZc/lb8YyDGX0cxoTvP9CMrj5dNsNmu/3YULY1Q2Tgx788+H2UwqEGyqMvNRXER8V5xRXBZvc8xh37DF9zn2E1vMsSfyF5THzLZJtScYo3FwEJks4Ky+ygUtrXLidAvsnA04wwzNKV2OLSurYXvFfbxAReJfgRIlew8PdGpF3qiDrvAnedLY2JOAJDYh8JNkuFDcb1naBFsgH6O68qYq6Nw49CJ3LbhNh5gGRCeBXE/Wfm+ZqNK2+SijQDeQc4WcRjLMq74VJWCJkp0M9aHv6trSBKkflJYGsfpyaDmqD/oBJcsQV9Vd23QTtSdWS37TEebgiJY1ruvBFQc2vam8Agnu4lfniKdlpxfvBy1fqlwJPQXQ2WpEkU03pZbGsaSCcBeRau2CmPcAtfUAz7S7rMOrVNy4+52n5lNMn3H7YhjDlkhSQldP15osszy8jQdDuPXJrsSzUuFdnU4ydF2u8IgRSMiu++fnpweDd4fFFu+Xjw+On+9y7DJofn0Y6IevRTb+a62mg/9tHtBNmMEz8oMMtl/ukfezFv4FhyamkA=="*)(*]VB*)}}

Plot it

Plot[Evaluate[y[x] /. s], {x, 0, 30}, PlotRange -> All]

(*VB[*)(FrontEndRef["3e8dfc8b-e98f-4890-97be-2e71c65f4da6"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKG6dapKQlWyTpplpapOmaWFga6FqaJ6XqGqWaGyabmaaZpCSaAQCQvBYk"*)(*]VB*)

Second-order nonlinear ordinary differential equation:

s = NDSolve[{y''[x] + Sin[y[x]] y[x] == 0, y[0] == 1, y'[0] == 0}, 
  y, {x, 0, 30}]

{{y->(*VB[*)BoxForm`temporalStorage$507286(*,*)(*"1:eJxdkUuTmkAAhM2rKpVfkezJha1SZPGRmzzkISA4MKJbe0AYcXgM6DCy+pfyJ6M55vJVd1efun/t6/XhS6/Xoz/u0FLc1meIUQf+fOr1TNKic1OXcYtJtmAkaXFN3voclN+455997oV77nNPwm9kXcswFOwPde5H+0E3HxqdSY6AOukyQMPNxa5GCm8pDYy8oGGW7DYFSC1hNptOaGzjONS2RuOItlherkbVgokm6bVItUTPh/7FKtbpzZpsrltKvSArKpQT2BA2c/19qA+Wp8mt7pZje3Fsk0QbXl2IpchbzPDwxurNuGN4tZ3qYbIWKkPNzCgKfOuYJyXUVcef57gyT2NDdqnmjEiaqGdpd2jLlOc7JBmvERPXtuukwIA395R5imEyEuyQYwaxBbAWq1GXpRlBikT0II63fCQqrWnByN/ly9VcTKBajHahJxlLku8/chwqtgAmfAUOgTk4y0AYu/wg9V5NxJ8UZ+TvYscFVJuKjJd5WWDmoY2mdN6hBSqUnHojfuobCayt4umx/zuUuef3w+fHh1/vWLMSge8PgeJ0RcrrvzQ4M/Rf53E4QCVK2nhfIvrtbhdxSdFfwbujog=="*)(*]VB*)}}

Plot[Evaluate[{y[x], y'[x], y''[x]} /. s], {x, 0, 30}, 
 PlotStyle -> Automatic, PlotLegends->Automatic]

(*VB[*)(Legended[ToExpression[FrontEndRef["0c3f6839-d6e3-4d6b-85c3-8d4d29103162"], InputForm], Placed[LineLegend[{Directive[Opacity[1.], RGBColor[0.24, 0.6, 0.8], AbsoluteThickness[2]], Directive[Opacity[1.], RGBColor[0.95, 0.627, 0.1425], AbsoluteThickness[2]], Directive[Opacity[1.], RGBColor[0.455, 0.7, 0.21], AbsoluteThickness[2]]}, {HoldForm[Placeholder[Style[1, Smaller]]], HoldForm[Placeholder[Style[2, Smaller]]], HoldForm[Placeholder[Style[3, Smaller]]]}, LegendMarkers -> None, LabelStyle -> {}, LegendLayout -> "Column"], After, Identity]])(*,*)(*"1:eJy1Uz1PG0EQPYOTAAlSiFIbIdFaAg4sU6ATxOZDssPHofR7t3OwYr2L5vYA0/MjoKWjoKampKCgoaKwEkVCSAj+AfuBLZkKkTDF093buTdvduZGIrmW9Hqelw5q+MVgtwKxRKIkhp80U4MNEHQi6TEpAxqqlOkzk5jkDPdNwzxKoaqCVvcgzhSJOISjmh6L/aRU9qeLtAR+cZKWomJ5KvaLZTpJJ6bHx/zx0rNwXsNapj/rMw9A6LLgTcuuYwbO30cNK5zEQJMPbTM1JsA5dDl5y6XKvfVrqDCEWLEdcG5NS8vbJGaqiZ6Nh8Al28oLcz8kl4hnhYO71bPLAH0bvwM8OjRxGziZIQ2zUSp5pmB9k8VbAtKUmVb+tXJi4z7Ak8erevT1T4AzA9fH2zPn71552LRcuAmeLfwNsLX//XShdfHKyp27z7WlFyWn8xIbjvncnt+mpgHd4M0kQ9XkwDpWwwbhXJ+/Tabn/8j0vpTpWlP7r7i9qxPcAkzt0U8p4EWiXVISAbeyidd1Ud2pXzqaNdKUmQrNwuvRZA1hjc0mShcynSxREEpP8wn5XfGE"*)(*]VB*)

The Lorenz equations

s = NDSolve[{Derivative[1][x][t] == -3 (x[t] - y[t]), 
   Derivative[1][y][t] == -x[t] z[t] + 26.5 x[t] - y[t], 
   Derivative[1][z][t] == x[t] y[t] - z[t], x[0] == z[0] == 0, 
   y[0] == 1}, {x, y, z}, {t, 0, 200}];

Plot it

ParametricPlot3D[Evaluate[{x[t], y[t], z[t]} /. s], {t, 0, 200}, 
 PlotPoints -> 1000]

(*VB[*)(FrontEndRef["91e247d8-8b60-4a10-8126-fa5c053d2056"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKWxqmGpmYp1joWiSZGeiaJBoa6FoYGpnppiWaJhuYGqcYGZiaAQB3JhTa"*)(*]VB*)

Simple wave evolution with periodic boundary conditions:

ClearAll[u]

s = NDSolve[{
  D[u[t,x], t,t] == D[u[t,x], x,x],
  u[0, x] == Exp[-(*SpB[*)Power[x(*|*),(*|*)2](*]SpB*)],
  Derivative[1,0][u][0,x] == 0, u[t, -10] == u[t, 10]
}, 
  u, {t, 0, 40}, {x, -10, 10}]

{{u->(*VB[*)BoxForm`temporalStorage$520073(*,*)(*"1:eJxdksmymlAARM1UlcpXJG/lI6lCEFCzA2WeFPCivHoLkMssF7jgwK/l56JZZtPV3dW70z8i5CSfJpMJ/vYQMc571IEcXt0/HyYTte5h16Aq7PM6lYb61OeofpsSQHgjXr9PiV/E65R4oX5D7V7t95Rx2/C7DpBM6gqM6F0bWrYrmppHnZcc4lI6tYbX0n45Vns9MsRzO5AkzGMJbNczlfGUcZ57o1KUFwewFxtxRZGmsFGyTjY03GUXAURUlUn2OgeZmcEbjNqdSKIWcwa2Wc1AEG9NRbDMBvIc697iowU6XcLYQKlNeu0g8vb2JzjgwuHXEXVbcxvRzAqFAfcrrFaRYGEZ7+OG1hakn7T0YrZKVxaZwYUm+znPhElNnZvrUt33AVJPNG9FwCgtQOmMEBwbMtWDZeXzZ/kW8JzFdrulS8YzBW5oVuwRclrUzYV+FI5roNyZ8XDRQ/1W5va1MBLDIwNn1WUwbcfDklqQ6lLOXIcXW6T6mnDofMlcsO2QgiSr6eLup2pozrltSevskPB8sFOj08uTxzsQiNf35OOT6eeHOEMF3a9PA8PYrqv7v9brBvjf5nkAF1bw1IdRBfGXR5TCCsO/oHyqMA=="*)(*]VB*)}}

Plot the solution:

DensityPlot[Evaluate[First[u[t, x] /. s]], {t, 0, 40}, {x, -10, 10}, 
 PlotPoints -> 30]

(*VB[*)(FrontEndRef["6e18cb89-fe5b-4cd2-b1ff-0c6df326f2fe"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKm6UaWiQnWVjqpqWaJumaJKcY6SYZpqXpGiSbpaQZG5mlGaWlAgCT/Rac"*)(*]VB*)